/*
  labFisObblStrut.cc

  codice C++ per Root di esempio per il lavoratorio di fisica SNS

  Il programma fa una simulazione Monte Carlo del rendimento
  di un'obbligazione strutturata delle Poste Italiane
  (http://www.poste.it/bancoposta/premium/scheda_premium.shtml).

  Si suppone che dato un il prezzo di un indice come f. del tempo X(t),
  ln[ X(t)/X(0) ] e' una variabile casuale Gaussiana con
  media mu*t/tau e varianza sigma^2*t/tau, dove

  - mu e' 1+r dove r e' l'interesse nel periodo tau
    - esempio 4% = 0.04 per un anno di interesse dei BOT

  - sigma e' la volatilita'
    - per un'azione puo' essere 1% al giorno

  Si suppone che tutti gli indici abbiano media che aumenta con
  uno stesso interesse, quello senza rischio dei titoli di Stato
  di breve durata, tipicamente 2%-4%.

  In questa obbligazione, la "cedola variabile" e' data da una certa
  percentuale [62%] della somma pesata dei tre indici (con peso
  maggiore [60%] per l'indice piu' alto, minore [40%] per il secondo e
  nullo per il piu' basso), e al cliente viene corrisposto il massimo
  tra la cedola variabile e un minimo garantito [8%].
  
  Questo programma produce un file labFisObblStrut2.root contenente
  tutti gli istogrammi generati.

  Come usarlo in ROOT (http://root.cern.ch):

    root labFisObblStrut2.cc+ <return>

  oppure

    root <return>
    root [0] .x labFisObblStrut.cc+

///////////////////////////////////////////////////////////////////////////////

  Revisioni:

  A.Lusiani, 15/5/2009
  - compila su Root 5.18

  A.Lusiani, 22/5/2009
  - heavily redesigned

*/

#if !defined(__CINT__) || defined(__MAKECINT__)

#include <assert.h>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <numeric>

#include <string.h>

#include "TROOT.h"
#include "TSystem.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TH1F.h"
#include "TMath.h"
#include "TFile.h"
#include "TTree.h"
#include "TRandom.h"

#endif

const char* rootFileName = "labFisObblStrut2.root";

void labFisObblStrut()
{
  gROOT->Reset();
  gROOT->SetStyle("Plain");

  //--- Create a new canvas.
  TCanvas* c1 = new TCanvas("c1", " ", 200, 10 , 700, 500);

  //
  // Create a new ROOT binary machine independent file.
  // This file is now becoming the current directory.
  //
  TFile* hfile =
    new TFile(rootFileName,"RECREATE","file obbligazione strutturata");
  
  //
  // create histograms for data accumulation
  //
  TH1F* ha = new TH1F("ha","Index A",100,-100,500);
  TH1F* hb = new TH1F("hb","Index B",100,-100,600);
  TH1F* hc = new TH1F("hc","Index C",100,-200,800);
  TH1F* hcedola = new TH1F("cedola","Cedola variabile",100,-100,500);
  TH1F* hfinal = new TH1F("final","Calcolo finale",100,-20,180);
  
  hcedola->Draw("hist");

  int nIter(10000);

  cout
    << "Simulazione Monte Carlo rendimento obbligazione strutturata"
    << ", " << nIter << " iterazioni"
    << endl;

  //
  // Parameters
  //

  //--- durata su quale si calcola il rendimento degli undici
  const int duration = 5;
  //--- numero di giorni all'anno (in prima approssimazione)
  const int daysPerYear = 365;

  //
  // tasso di interesse annuale senza rischio in per cento
  // (uguale per tutti gli indici)
  //
  const float interest = 4;
  //--- tasso di accrescimento annuale senza rischio
  const float rate = log(1+interest/100);

  //--- volatilita' giornaliera dei tre indici
  const float vol_a = float(0.5)/100;
  const float vol_b = float(1)/100;
  const float vol_c = float(2)/100;

  //--- minimo guadagno garantito in per cento
  const float garantito = 8;
  //--- quota di partecipazione al guadagno degli indici
  const float partecipazione = 0.62;
  
  const int kUPDATE = 100;

  //
  // Variables
  //

  double a, b, c;
  double sum, cedola_variabile, final;
  double meanGain(0);
  double meanGainErr(0);

  //--- the seed determines the pseudo-random sequence
  gRandom->SetSeed(1);

  // Loop over "experiments":
  for (int i=0; i < nIter; i++) {

    // --- extract ln[ X(duration) / X(0) ] for the 3 indices
    a = gRandom->Gaus( duration*rate , sqrt(duration*daysPerYear)*vol_a );
    b = gRandom->Gaus( duration*rate , sqrt(duration*daysPerYear)*vol_b );
    c = gRandom->Gaus( duration*rate , sqrt(duration*daysPerYear)*vol_c );

    // cout << a << " " << b << " " << c << endl;

    //--- invert the logarithm to get  X(duration) / X(0) 
    a = exp(a);
    b = exp(b);
    c = exp(c);

    //--- get the percent gain
    a = (a-1)*100;
    b = (b-1)*100;
    c = (c-1)*100;
    
    //--- accumulate the distribution of percent gains
    ha->Fill(a);
    hb->Fill(b);
    hc->Fill(c);

    //--- sort percent gains in ascending order
    vector<double> sortedgains(3);
    sortedgains[0] = a;
    sortedgains[1] = b;
    sortedgains[2] = c;
    sort( sortedgains.begin(), sortedgains.end() );

#if 0
    cout
      << sortedgains[0] << " "
      << sortedgains[1] << " "
      << sortedgains[2] << endl;
#endif

    //--- compute weighed sum of gains and the client gain
    sum = 0.6*sortedgains[2] + 0.4*sortedgains[1];
    cedola_variabile = partecipazione*sum;
    hcedola->Fill(cedola_variabile);

    //--- give at least minimum gain
    final = max(double(garantito),cedola_variabile);
    hfinal->Fill(final);

    //--- accumulate gains and gains^2
    meanGain += final;
    meanGainErr += pow(final,2);

    if ((i%kUPDATE) == (kUPDATE-1)) {
      //--- update display
      c1->Modified();
      c1->Update();
    }
  }

  //--- update display
  c1->Modified();
  c1->Update();
  
  //--- save all objects in this file
  hfile->Write();
  
  //--- print average gain in percent and its error
  meanGain = meanGain/nIter;
  meanGainErr = meanGainErr - nIter*pow(meanGain,2);
  meanGainErr = sqrt(meanGainErr)/nIter;
  cout
    << "mean gain: " << meanGain << " +- " << meanGainErr << endl;

  //
  // print fair price
  // exp(-rate*t/tau) * price(t) = 
  // = exp(-rate*t/tau) price(0) * (1+percent_gain/100)
  //
  cout
    << "fair price: " << exp(-rate*duration) * (1+meanGain/100) << endl;

  //
  // Note that the file is automatically closed when application terminates
  // or when the file destructor is called.
  //
}
